C     MINSUB      Minimization of the sum of squares.
C     Copyright (C) 1997  Filip Hroch, Masaryk University, Brno, CZ
C                    
C     This program is free software; you can redistribute it and/or modify
C     it under the terms of the GNU General Public License as published by
C     the Free Software Foundation; either version 2 of the License, or
C     (at your option) any later version.
C     
C     This program is distributed in the hope that it will be useful,
C     but WITHOUT ANY WARRANTY; without even the implied warranty of
C     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C     GNU General Public License for more details.
C
C     You should have received a copy of the GNU General Public License  
C     along with this program; if not, write to the Free Software
C     Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
C
C    E-mail: hroch@physics.muni.cz
C    Web:    http://www.sci.muni.cz/~hroch/index.html
C    Postal address: F.Hroch
C                    Department of Theoretical Physics and Astrophysics
C                    Masaryk University
C                    Kotlarska 2
C                    611 37 Brno
C                    Czech Republic
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      FUNCTION TRAFO(NS,CLIP,ID1,ID2,N1,N2,X1,Y1,X2,Y2,R,M,C,N)
C Find coefficients of transformation.
C 
      DOUBLE PRECISION TRAFO, MACHMAX,MACHEPS,ARC1
      PARAMETER( MACHMAX = 1D10, MACHEPS = 1D-6, ARC1 = 0.2 )

      INTEGER ID1(*),ID2(*),I,RESROV,J,J0,N,NS,NN,N1,N2
      DOUBLE PRECISION R(2),M(2,2),X1(*),Y1(*),X2(*),Y2(*),C,CLIP
      DOUBLE PRECISION B(2),A(2),MIROR,S,XX1,YY1,XX2,YY2,XX,YY
      DOUBLE PRECISION DR,DX,DY,TOLS,RMIN
      DOUBLE PRECISION AXI,AYI,A2,P(4),D(4,4),X(4)

      TRAFO = MACHMAX
C Mirroring
      A(1) = X1(ID1(2))-X1(ID1(1))
      A(2) = Y1(ID1(2))-Y1(ID1(1))
      B(1) = X1(ID1(3))-X1(ID1(1))
      B(2) = Y1(ID1(3))-Y1(ID1(1))
      MIROR = A(1)*B(2)-A(2)*B(1)
      A(1) = X2(ID2(2))-X2(ID2(1))
      A(2) = Y2(ID2(2))-Y2(ID2(1))
      B(1) = X2(ID2(3))-X2(ID2(1))
      B(2) = Y2(ID2(3))-Y2(ID2(1))
      MIROR = MIROR*(A(1)*B(2)-A(2)*B(1))
      MIROR = MIROR/ABS(MIROR)
C Zero-iteration
      N = 0
      DO I=1,4
        P(I) = 0D0
        DO J=1,4
          D(I,J) = 0D0
        ENDDO
      ENDDO
      AXI = 0D0
      AYI = 0D0
      A2 = 0D0
      DO I=1,NS
        XX1 = X1(ID1(I))
        YY1 = Y1(ID1(I))
        XX2 = X2(ID2(I))
        YY2 = Y2(ID2(I))
        AXI = AXI+XX1
        AYI = AYI+YY1
        A2 = A2+XX1**2+YY1**2
        P(1) = P(1)+XX2
        P(2) = P(2)+YY2
        P(3) = P(3)+(XX1*XX2+MIROR*YY1*YY2)
        P(4) = P(4)+(YY1*XX2-MIROR*XX1*YY2)
      ENDDO
      D(1,1) = NS*C
      D(1,3) = AXI*C
      D(1,4) = AYI*C
      D(2,2) = NS*C
      D(2,3) = MIROR*AYI*C
      D(2,4) = -MIROR*AXI*C
      D(3,1) = AXI*C
      D(3,2) = MIROR*AYI*C
      D(3,3) = A2*C
      D(4,1) = AYI*C
      D(4,2) = -MIROR*AXI*C
      D(4,4) = D(3,3)
      IF( RESROV(4,4,D,P,X).NE.0 ) RETURN
C Goniometric unit.
c      print *,'1>',X(3)**2+X(4)**2
      IF( ABS(X(3)**2+X(4)**2-1D0).GT.ARC1) RETURN
      R(1) = X(1)
      R(2) = X(2)
      M(1,1) = X(3)
      M(1,2) = X(4)
      M(2,1) = -MIROR*X(4)
      M(2,2) = MIROR*X(3)
C Reziduals.
      S = 0D0
      DO I=1,NS
        XX1 = X1(ID1(I))
        YY1 = Y1(ID1(I))
        XX = C*(R(1)+XX1*M(1,1)+YY1*M(1,2))
        YY = C*(R(2)+XX1*M(2,1)+YY1*M(2,2))
        XX2 = X2(ID2(I))
        YY2 = Y2(ID2(I))
C        print *,i,XX-XX2,YY-YY2
        S = S+(XX-XX2)**2+(YY-YY2)**2
      ENDDO
      TOLS = CLIP**2*S/(NS-4D0)+MACHEPS
C      print *,S,sqrt(TOLS),c

C Rezidulas for all stars
      NN = 0
      S = 0D0
      DO I=1,4
        P(I) = 0D0
        DO J=1,4
          D(I,J) = 0D0
        ENDDO
      ENDDO
      AXI = 0D0
      AYI = 0D0
      A2 = 0D0
C  
      DO I=1,N1
        XX = C*(R(1)+X1(I)*M(1,1)+Y1(I)*M(1,2))
        YY = C*(R(2)+X1(I)*M(2,1)+Y1(I)*M(2,2))
C Find star with the smallest residues
        RMIN = MACHMAX
        DO J=1,N2
          DX = XX-X2(J)
          DY = YY-Y2(J)
          DR = DX**2+DY**2
          IF( DR.LT.RMIN )THEN
            RMIN = DR
            J0 = J
          ENDIF
        ENDDO
C 
        J = J0
c        print *,'.',i,j,rmin,tols
        IF( RMIN.LE.TOLS )THEN
          S = S+RMIN
          NN = NN+1
          AXI = AXI+X1(I)
          AYI = AYI+Y1(I)
          A2 = A2+X1(I)**2+Y1(I)**2
          P(1) = P(1)+X2(J)
          P(2) = P(2)+Y2(J)
          P(3) = P(3)+(X1(I)*X2(J)+MIROR*Y1(I)*Y2(J))
          P(4) = P(4)+(Y1(I)*X2(J)-MIROR*X1(I)*Y2(J))
c           PRINT *,I,J,XX-X2(J),YY-Y2(J)
        ENDIF
      ENDDO
      IF(NN.LE.4) RETURN
      TOLS = CLIP**2*S/(NN-4D0)+MACHEPS
      D(1,1) = NN*C
      D(1,3) = C*AXI
      D(1,4) = C*AYI
      D(2,2) = NN*C
      D(2,3) = MIROR*C*AYI
      D(2,4) = -MIROR*C*AXI
      D(3,1) = C*AXI
      D(3,2) = MIROR*C*AYI
      D(3,3) = C*A2
      D(4,1) = C*AYI
      D(4,2) = -MIROR*C*AXI
      D(4,4) = D(3,3)
      IF( RESROV(4,4,D,P,X).NE.0 ) RETURN
C  Goniometric unit.
c      print *,'2>',X(3)**2+X(4)**2
      IF( ABS(X(3)**2+X(4)**2-1D0).GT.ARC1) RETURN
C      PRINT *,'>',NN,(X(I),I=1,4),TOLS
      R(1) = X(1)
      R(2) = X(2)
      M(1,1) = X(3)
      M(1,2) = X(4)
      M(2,1) = -MIROR*X(4)
      M(2,2) = MIROR*X(3)

      NN = 0
      S = 0D0
C TEST SPRAVNOSTI VYSLEDNE ZPRESNENE TRANSORMACE VYPOCET VYSLEDNE SUMY CTVERCU
C A POCTU HVEZD
      DO I=1,N1
        XX = C*(R(1)+X1(I)*M(1,1)+Y1(I)*M(1,2))
        YY = C*(R(2)+X1(I)*M(2,1)+Y1(I)*M(2,2))
        RMIN = MACHMAX
        DO J=1,N2
          DX = XX-X2(J)
          DY = YY-Y2(J)
          RMIN = MIN(RMIN,DX**2+DY**2)
        ENDDO
c        print *,i,j,rmin,tols
        IF( RMIN.LT.TOLS )THEN
          S = S+RMIN
          NN = NN+1
        ENDIF
      ENDDO
      IF( NN.GT.4 ) THEN
        TRAFO = S/(NN - 4.0)+MACHEPS
        N = NN
      ENDIF

      END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      INTEGER FUNCTION RESROV(NMAX,N,A,B,X)

      DOUBLE PRECISION A(4,4),B(4),X(4),T
      DOUBLE PRECISION A1(5,4)
      INTEGER NMAX,N,I,J

      RESROV = 1
      DO I=1,N
        DO J=I,N
          A1(I,J) = A(I,J)
        ENDDO
      ENDDO
      CALL DLL(4,4,A1,T)
      IF( ABS(T).LT.1E-5 ) RETURN
      CALL DLLSOL(4,4,A1,B,X)
      RESROV = 0
      END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
